SoilBalance.f90 Source File
Solve soil water and energy balance
Source Code
!! Solve soil water and energy balance !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.6 - 20th July 2023 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 28/Jun/2010 | Initial code. chiara corbari and davide rabuffetti | ! | 1.1 | 03/May/2011 | Correct energy fluxes when snow is present | ! | 1.2 | 10/Nov/2016 | Module rewritten to add new infiltration models and different soil properties management | ! | 1.3 | 25/Jul/2022 | subsurface flow changed from muskingum to darcy. Runon option removed | ! | 1.4 | 27/Jan/2023 | bug correction: percolationcell set to zero on hillslope and rivers in SetSoilDepth | ! | 1.5 | 10/Feb/2023 | modified to not use river network | ! | 1.6 | 20/Jul/2023 | SetSoilDepth renamed to PercolationAndCaprise and modified | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Solve soil water and energy balance to update soil moisture and ! compute actual evapotranspiration, runoff, and groundwater ! recharge. ! Soil column is assumed to be divided in two layers: ! ! * the surface layer where roots develop (root zone) ! ! * the lower unsaturated layer that transfers ! water to recharge groundwater in landplain (transmission zone) ! MODULE SoilBalance ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Type Definitions: short, float, double USE DomainProperties, ONLY : & !imported variables: mask USE StringManipulation, ONLY: & !Imported routines: ToString USE GridLib, ONLY: & !Imported definitions: grid_integer, grid_real, & !Imported routines: NewGrid USE GridOperations, ONLY : & !imported routines: GridByIni, CRSisEqual, CellArea, & !Imported assignment: ASSIGNMENT (=) USE Chronos, ONLY: & !Imported definitions: DateTime, & !Imported operators and assignment: OPERATOR (>=), OPERATOR (<=), & ASSIGNMENT( = ) USE LogLib, ONLY: & !Imported routines: Catch USE IniLib, ONLY: & !Imported derived types: IniList, & !Imported routines: IniOpen, IniClose, & IniReadInt, IniReadReal, IniReadDouble, & IniReadString, SectionIsPresent, & KeyIsPresent USE SoilProperties, ONLY : & !Imported definitions: !SoilType, & !imported routines: UnsHydCond, Psi USE Infiltration, ONLY : & !imported routines: InfiltrationInit, & Sat2scn, SCScurveNumber, & Philip, GreenAmpt, & !imported variables: infiltrationModel, & soilDivisions, & curveNumber, & storativity, & sEff, ism, cuminf, & raincum, wiltingPoint, & fieldCapacity, thetar, & thetas, ksat, psdi, psic, & abstractionRatio, phy, & !imported parameters: SCS_CN, PHILIPEQ, & GREEN_AMPT, ROSS_BC, ROSS_VG USE Richards, ONLY : & !Imported routines SolveRichardsBC, & !imported variables : nsteps, parameters, cell, & dx, S, h0grid, drncum, & infilcum, evapcum, runoffcum USE Evapotranspiration, ONLY : & !imported routines ETinit, PETupdate, & !imported variables: pet, pe, pt USE Morphology, ONLY: & !Imported routines: DownstreamCell, CheckOutlet USE TableLib, ONLY: & ! Imported definitions: Table, & !Imported routines TableNew, TableGetValue, TableGetNrows USE StringManipulation, ONLY: & !Imported routines: StringToShort USE MorphologicalProperties, ONLY: & !imported variables: dem USE Snow, ONLY : & !imported variables: swe, waterInSnow, dtSnow USE Glacier, ONLY : & !imported variables: icewe, waterInIce, dtIce !Global declarations: !public routines PUBLIC :: InitSoilBalance PUBLIC :: SolveSoilBalance !public variables: INTEGER (KIND = short) :: dtSoilBalance !TYPE(NetworkSoilBalance):: soilNetwork !!define spatial scheme to solve balance !parameters TYPE (grid_real) :: ksat_sub ![m/s] TYPE (grid_real) :: soilDepth !!soil depth [m] TYPE (grid_real) :: soilDepthRZ !! root zone depth [m] TYPE (grid_real) :: soilDepthTZ !!transmission zone depth [m] !common state variables TYPE (grid_real) :: soilSat !!mean soil relative saturation [0-1] TYPE (grid_real) :: soilSatRZ !!root zone soil relative saturation [0-1] TYPE (grid_real) :: soilSatTZ !!transmission zone soil relative saturation [0-1] TYPE (grid_real) :: soilMoisture !!mean volumetric water content [m3/m3] TYPE (grid_real) :: soilMoistureRZ !!root zone volumetric water content [m3/m3] TYPE (grid_real) :: soilMoistureTZ !!transmission zone volumetric water content [m3/m3] TYPE (grid_integer) :: rainFlag !!set if it is raining in a cell TYPE (grid_integer) :: interstormDuration !!interstorm duration [s] TYPE (grid_real) :: QinSoilSub !!input discharge in soil subsurface [m3/s] TYPE (grid_real) :: QoutSoilSub !!output discharge from soil subsurface [m3/s] TYPE (grid_real) :: balanceError !!balance error [mm] TYPE (grid_real) :: deltaSoilMoisture !! time step soil mositure variation (m) !vertical fluxes TYPE (grid_real) :: rainBalance !!actual rainfall rate as input to soil balance [m/s] TYPE (grid_real) :: infilt !!infiltration rate [m/s] TYPE (grid_real) :: runoff !!runoff rate [m/s] TYPE (grid_real) :: et !! actual evapotranspiration rate [m/s] TYPE (grid_real) :: percolation !! deep percolation through soil [m/s] TYPE (grid_real) :: percolationFactor !! deep percolation factor [-] TYPE (grid_real) :: capRise !! capillary rise [m/s] !energy balance fluxes TYPE (grid_real) :: Ts TYPE (grid_real) :: Rnetta TYPE (grid_real) :: Xle TYPE (grid_real) :: Hse TYPE (grid_real) :: Ge TYPE (grid_real) :: Ta_prec TYPE (grid_real) :: Ts_prec !Local routines PRIVATE :: WetlandIsWet PRIVATE :: SetInitialCondition PRIVATE :: SetWetland PRIVATE :: PercolationAndCaprise PRIVATE :: UpdateSoilMoisture !local declarations: !parameters INTEGER, PRIVATE, PARAMETER :: HILLSLOPE = 0 INTEGER, PRIVATE, PARAMETER :: CHANNEL = 1 INTEGER, PRIVATE, PARAMETER :: LAKE = 2 INTEGER, PRIVATE, PARAMETER :: LANDPLAIN = 3 !wetland TYPE (grid_integer), PRIVATE :: wetland INTEGER, ALLOCATABLE, PRIVATE :: wetCode (:) TYPE (DateTime), ALLOCATABLE, PRIVATE :: wetBegin (:) TYPE (DateTime), ALLOCATABLE, PRIVATE :: wetEnd (:) !balance REAL (KIND = float), PRIVATE :: thresholdStartEvent !threshold to consider storm initiated [m/s] REAL (KIND = float), PRIVATE :: interstorm !duration of interstorm period to terminate an event [s] REAL (KIND = float), PRIVATE :: isd !!initial saturation degree, used for cold start TYPE (grid_integer), PRIVATE :: balanceId !!id code for solving water balance LOGICAL , PRIVATE :: saturatedByGroundwater !!groundwater table intercepts root zone !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize soil water balance SUBROUTINE InitSoilBalance & ! (inifile, flowDirection, time) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information TYPE (grid_integer), INTENT(IN) :: flowDirection !!flow direction map TYPE (DateTime), INTENT(IN) :: time !!start time !local declarations: TYPE (IniList) :: iniDB CHARACTER (LEN = 1000) :: filename INTEGER (KIND = short) :: k, i, j, iin, jin, is, js REAL (KIND = float) :: scalar !------------end of declaration------------------------------------------------ !open and read configuration file CALL IniOpen (inifile, iniDB) !read balance id IF (SectionIsPresent('balance-id', iniDB)) THEN CALL GridByIni (iniDB, balanceId, section = 'balance-id') ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading balance-id: ' , & argument = 'section not defined in ini file' ) END IF !read subsurface soil hydraulic conductivity IF (SectionIsPresent('ksat-subsurface', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'ksat-subsurface') ) THEN scalar = IniReadReal ('scalar', iniDB, 'ksat-subsurface') CALL NewGrid (ksat_sub, mask, scalar) ELSE CALL GridByIni (iniDB, ksat_sub, section = 'ksat-subsurface') END IF ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading subsurface conductivity: ' , & argument = 'section not defined in ini file' ) END IF !read soil depth (m) IF (SectionIsPresent('soil-depth', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'soil-depth') ) THEN scalar = IniReadReal ('scalar', iniDB, 'soil-depth') CALL NewGrid (soilDepth, mask, scalar) ELSE CALL GridByIni (iniDB, soilDepth, section = 'soil-depth') END IF ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading soil depth: ' , & argument = 'section not defined in ini file' ) END IF !read deep percolation factor (-) IF (SectionIsPresent('deep-percolation-factor', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'deep-percolation-factor') ) THEN scalar = IniReadReal ('scalar', iniDB, 'deep-percolation-factor') CALL NewGrid (percolationFactor, mask, scalar) ELSE CALL GridByIni (iniDB, percolationFactor, section = 'deep-percolation-factor') END IF ELSE !grid is optional: default = 1. CALL NewGrid (percolationFactor, mask, 1.0) END IF !read root zone depth (m) IF (SectionIsPresent('root-zone-depth', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'root-zone-depth') ) THEN scalar = IniReadReal ('scalar', iniDB, 'root-zone-depth') CALL NewGrid (soilDepthRZ, mask, scalar) ELSE CALL GridByIni (iniDB, soilDepthRZ, section = 'root-zone-depth') END IF ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading root zone depth: ' , & argument = 'section not defined in ini file' ) END IF !compute transmission zone depth CALL NewGrid ( soilDepthTZ, mask, 0.) DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) /= mask % nodata ) THEN soilDepthTZ % mat(i,j) = soilDepth % mat(i,j) - & soilDepthRZ % mat(i,j) IF ( soilDepthTZ % mat(i,j) <= 0. ) THEN soilDepthTZ % mat(i,j) = 0.1 END IF END IF END DO END DO !configure evapotranspiration IF (KeyIsPresent('evapotranspiration', iniDB)) THEN filename = IniReadString ('evapotranspiration', iniDB) CALL ETinit (filename, time) ELSE CALL Catch ('error', 'SoilBalance', & 'evapotranspiration not found in configuration file' ) END IF !threshold to initiate storm period, read and convert from mm/h to m/s thresholdStartEvent = IniReadReal ('threshold-storm-start', iniDB) thresholdStartEvent = thresholdStartEvent / 1000. / 3600. !interstorm duration to terminate an event, read and convert from hours to seconds interstorm = IniReadReal ('interstorm', iniDB) interstorm = interstorm * 3600. !configure infiltration IF (KeyIsPresent('infiltration', iniDB)) THEN filename = IniReadString ('infiltration', iniDB) CALL InfiltrationInit (filename, soilMoisture, soilDepth ) ELSE CALL Catch ('error', 'SoilBalance', & 'infiltration not found in configuration file' ) END IF !set initial condition CALL SetInitialCondition (iniDB) !set wetland CALL SetWetland (iniDB, inifile) !allocate coomon variables used to solve soil water and energy balance !vertical fluxes CALL NewGrid (rainBalance, mask, 0.) CALL NewGrid (infilt, mask, 0.) CALL NewGrid (runoff, mask, 0.) CALL NewGrid (et, mask, 0.) CALL NewGrid (percolation, mask, 0.) CALL NewGrid (capRise, mask, 0.) !energy balance CALL NewGrid (Ts, mask, 0.) CALL NewGrid (Rnetta, mask, 0.) CALL NewGrid (Xle, mask, 0.) CALL NewGrid (Hse, mask, 0.) CALL NewGrid (Ge, mask, 0.) CALL NewGrid (Ta_prec, mask, 0.) CALL NewGrid (Ts_prec, mask, 0.) !lateral fluxes at time t CALL NewGrid (QinSoilSub, mask, 0.) CALL NewGrid (QoutSoilSub, mask, 0.) !balance error CALL NewGrid (balanceError, mask, 0.) !soil mositure variation CALL NewGrid (deltaSoilMoisture, mask, 0.) ! Configuration terminated. Deallocate ini database CALL IniClose (iniDB) RETURN END SUBROUTINE InitSoilBalance !============================================================================== !| Description: ! Solve soil water balance SUBROUTINE SolveSoilBalance & ! (time, rain, irrigation, flowdir, vf, vadose) IMPLICIT NONE !Arguments with intent(in): !INTEGER (KIND = short), INTENT(IN) :: dt !!computation time step [s] TYPE (DateTime), INTENT(IN) :: time !current date and time TYPE (grid_real), INTENT(IN) :: rain !!rainfall rate [m/s] TYPE (grid_real), INTENT(IN) :: irrigation !!irrigation rate [m/s] TYPE (grid_integer), INTENT(IN) :: flowdir TYPE (grid_real), INTENT(IN) :: vf !!vegetation fraction [0-1] TYPE (grid_real), INTENT(IN) :: vadose !vadose zone depth [m] !local declarations: INTEGER (KIND = short) :: k, iin, jin, is, js, id, i, j REAL (KIND = double) :: runoffcell !runoff of single cell [m/s] REAL (KIND = float) :: percolationcellRZ !root zone percolation of single cell [m/s] REAL (KIND = float) :: percolationcellTZ !transmission zone percolation of single cell [m/s] REAL (KIND = float) :: Tdeep !deep surface temperature in lake REAL (KIND = float) :: alpha, beta !!used to compute actual evapotranspiration [0-1] REAL (KIND = float) :: sm, wp, fc !!soil moisture, wilting point and field capacity of current cell REAL (KIND = float) :: actualE, actualT !!actual evaporation and transpiration of current cell REAL (KIND = double) :: soilDepthCell !!soil depth of current cell [m] REAL (KIND = float) :: prunoffcum !!cumulated runoff computed by Ross at previous timestep [cm] REAL (KIND = float) :: pinfilcum !!cumulated infiltration computed by Ross at previous timestep [cm] REAL (KIND = float) :: pdrncum !!cumulated drainage computed by Ross at previous timestep [cm] REAL (KIND = float) :: pevapcum !!cumulated evapotranspiration used and possibly modified by Ross at previous timestep [cm] REAL (KIND = float) :: water !!water in soil [m] REAL (KIND = float) :: ddx !actual cell length (m) REAL (KIND = float) :: meanHydCond !!mean hydraulic conductivity !for computing lateral subsurface flow (m/s) REAL (KIND = float) :: ds !!thickness of the saturated depth (m) REAL (KIND = float) :: cellWidth !!cell width (m) REAL (KIND = float) :: ijSlope LOGICAL :: snowCovered, iceCovered !------------end of declaration------------------------------------------------ !reset variables for the time step QinSoilSub = 0. QoutSoilSub = 0. runoff = 0. rainBalance = 0. !update PET !---------- CALL PETupdate (time) !compute lateral subsurface intercell fluxes !------------------------------------------ DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) == mask % nodata ) CYCLE !check balance id IF ( balanceId % mat (i,j) == LANDPLAIN .OR. & balanceId % mat (i,j) == LAKE ) THEN CYCLE END IF !set downstream cell (is,js) CALL DownstreamCell(i, j, flowdir % mat(i,j), is, js, ddx, flowdir ) IF ( CheckOutlet (i, j, is, js, flowdir) ) CYCLE !harmonic mean saturated conductivity meanHydCond = 2. * ksat_sub % mat (i,j) * ksat_sub % mat (is,js) / & ( ksat_sub % mat (i,j) + ksat_sub % mat (is,js) ) !thickness of the saturated depth ds = soilDepthTZ % mat(i,j) * soilMoistureTZ % mat (i,j) / thetas % mat (i,j) !cell width cellWidth = ( CellArea (soilMoisture, i, j) ) ** 0.5 !local slope ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx IF ( ijSlope <= 0. ) THEN ijSlope = 0.0001 END IF QoutSoilSub % mat (i,j) = cellWidth * ds * meanHydCond * ijSlope !output becomes input of downstream cell QinSoilSub % mat (is,js) = QinSoilSub % mat (is,js) + QoutSoilSub % mat (i,j) END DO END DO !loop through mask DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) == mask % nodata ) CYCLE !reset fluxes of current cell runoffcell = 0. percolationcellRZ = 0. percolationcellTZ = 0. !set balance id id = balanceId % mat (i,j) !rain contributes to soil balance rainBalance % mat (i,j) = rain % mat (i,j) !add irrigation IF ( ALLOCATED ( irrigation % mat ) ) THEN rainBalance % mat (i,j) = rainBalance % mat (i,j) + irrigation % mat (i,j) END IF snowCovered = .FALSE. iceCovered = .FALSE. !check if snow and ice cover the current cell IF ( dtIce > 0 ) THEN !glaciers are simulated IF ( icewe % mat (i,j) <= 0. ) THEN !liquid water in ice !contributes to soil balance rainBalance % mat (i,j) = rainBalance % mat (i,j) + & waterInIce % mat (i,j) / dtSoilBalance waterInIce % mat (i,j) = 0. IF ( dtSnow > 0 ) THEN !snow is simulated IF ( swe % mat (i,j) <= 0. ) THEN !liquid water in snow ! contributes to soil balance rainBalance % mat (i,j) = rainBalance % mat (i,j) + & waterInSnow % mat (i,j) / dtSoilBalance waterInSnow % mat (i,j) = 0. ELSE rainBalance % mat (i,j) = 0. snowCovered = .TRUE. END IF END IF ELSE !soil is frozen, no vertical fluxes in it !runoff % mat(i,j) = 0. !infilt % mat(i,j) = 0. !percolation % mat(i,j) = 0. !et % mat(i,j) = 0. rainBalance % mat (i,j) = 0. iceCovered = .TRUE. END IF END IF IF ( dtSnow > 0 .AND. dtIce <= 0 ) THEN IF ( swe % mat (i,j) > 0.) THEN !soil is frozen, no vertical fluxes in it !runoff % mat(i,j) = 0. !infilt % mat(i,j) = 0. !percolation % mat(i,j) = 0. !et % mat(i,j) = 0. rainBalance % mat (i,j) = 0. snowCovered = .TRUE. ELSE rainBalance % mat (i,j) = rainBalance % mat (i,j) + & waterInSnow % mat (i,j) / dtSoilBalance waterInSnow % mat (i,j) = 0. END IF END IF !decide if storm or interstorm period IF ( rainFlag % mat(i,j) == 0 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN !storm starts rainFlag % mat(i,j) = 1 interstormDuration % mat(i,j) = 0 !reset interstorm duration IF (infiltrationModel == SCS_CN) THEN CALL Sat2scn (curveNumber % mat(i,j), storativity % mat(i,j), & soilSatRZ % mat(i,j), sEff % mat(i,j) ) END IF IF (infiltrationModel == PHILIPEQ .OR. & infiltrationModel == GREEN_AMPT) THEN ism % mat(i,j) = soilMoistureRZ % mat(i,j) END IF IF (infiltrationModel == ROSS_BC .OR. & infiltrationModel == ROSS_VG) THEN runoffcum % mat(i,j) = 0. drncum % mat(i,j) = 0. infilcum % mat(i,j) = 0. evapcum % mat(i,j) = 0. END IF ELSE IF ( rainFlag % mat(i,j) == 1 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN !storm continues rainFlag % mat(i,j) = 2 ELSE IF ( rainBalance % mat (i, j) < thresholdStartEvent ) THEN !update interstorm duration interstormDuration % mat(i,j) = interstormDuration % mat(i,j) + dtSoilBalance END IF !check if storm event must be terminated IF (interstormDuration % mat(i,j) >= interstorm) THEN rainFlag % mat(i,j) = 0 IF (infiltrationModel == SCS_CN) THEN !reset cumulated rainfall raincum % mat(i,j) = 0. END IF IF (infiltrationModel == PHILIPEQ .OR. & infiltrationModel == GREEN_AMPT ) THEN !reset cumulated infiltration cuminf % mat(i,j) = 0. END IF END IF !compute vertical fluxes according to soil balance computation scheme SELECT CASE ( id ) CASE (LAKE) !Lake drains subsurface flux !runoff can become < 0 when pet > rain + QinSoilSub runoffcell = rainBalance % mat (i, j) - pet % mat(i,j) + & QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j) et % mat(i,j) = pet % mat(i,j) !et is potential percolationcellRZ = 0. ! percolation is zero percolationcellTZ = 0. !CASE (LAKE_EWB) ! Tdeep = 25. ![°C] ! GR TODO !CALL ET_energy_nr_lago(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),& ! dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),& ! radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),& ! indice_nuv%mat(iin,jin),ET%mat(iin,jin),& ! Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),& ! Ge%mat(iin,jin),Tdeep) !runoffcell = rainBalance % mat (i, j) - et % mat(i,j) + & ! QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j) ! percolationcell = 0. !deep percolation is zero CASE DEFAULT ! all other cells(SLOPE, SLOPE_EWB, SLOPE_EWB_WET, CHANNEL, CHANNEL_EWB, & !CHANNEL_EWB_WET, PLAIN, PLAIN_EWB, PLAIN_EWB_WET) !compute infiltration and runoff during storm period IF ( rainFlag % mat(i,j) /= 0 ) THEN SELECT CASE (infiltrationModel) CASE (SCS_CN) !use soil conservation services curve number method infilt % mat (i,j) = SCScurveNumber (rain = rainBalance % mat (i, j), & raincum = raincum % mat (i,j), & c = abstractionRatio % mat(i,j), & s = sEff % mat(i,j), & dt = dtSoilBalance, runoff = runoffcell) CASE (PHILIPEQ) !use Philips infiltration model infilt % mat (i,j) = Philip (rain = rainBalance % mat (i, j), & psic = psic % mat(i,j), & ksat = ksat % mat(i,j), & theta = soilMoistureRZ % mat(i,j), & thetas = thetas % mat(i,j), & thetar = thetar % mat(i,j), & psdi = psdi % mat (i,j), & itheta = ism % mat (i,j), dt = dtSoilBalance, & cuminf = cuminf % mat (i,j) ) runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j) CASE (GREEN_AMPT) !use Green-Ampt infiltration model infilt % mat (i,j) = GreenAmpt (rainFlag % mat(i,j), & rain = rainBalance % mat (i, j), & ksat = ksat % mat(i,j), & theta = soilMoistureRZ % mat(i,j), & thetas = thetas % mat(i,j), & thetar = thetar % mat(i,j), & phy = phy % mat(i,j), & itheta = ism % mat (i,j), dt = dtSoilBalance, & cuminf = cuminf % mat (i,j) ) runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j) !CASE (ROSS_BC) !the subroutine SolveRichardsBC is always called even during interstorm period END SELECT ELSE !interstorm period, runoff is zero. Low intensity rainfall infiltrates runoffcell = 0. infilt % mat(i,j) = rainBalance % mat (i, j) END IF !compute actual evapotranspiration !--------------------------------- !IF (id == SLOPE_EWB .OR. id == SLOPE_EWB_WET .OR. id == CHANNEL_EWB .OR. & ! id == CHANNEL_EWB_WET .OR. id == PLAIN_EWB .OR. id == PLAIN_EWB_WET) THEN !solve energy balance to compute actual evapotranspiration !IF ( WetlandIsWet (wetland, iin, jin, time) ) THEN ! TODO ! Tsprof=26.5 ! !CALL ET_energy_nr_wetland(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),& ! dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),& ! radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),& ! indice_nuv%mat(iin,jin),ET%mat(iin,jin),& ! Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),& ! Ge%mat(iin,jin),Tsprof) ! ! if(ET%mat(iin,jin)<0)then ! ET%mat(iin,jin)=0 ! end if !ELSE ! CALL ET_ENERGETICO_suolo(SM%mat(iin,jin),SMsat%mat(iin,jin),SMres%mat(iin,jin),albedo%mat(iin,jin),& ! res_stom_min%mat(iin,jin),LAI%mat(iin,jin),h_veg%mat(iin,jin),bpress%mat(iin,jin),& ! dtm%mat(iin,jin),B%mat(iin,jin),fraz_veg%mat(iin,jin),temperatura_oss%mat(iin,jin),& ! radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),& ! indice_nuv%mat(iin,jin),WP%mat(iin,jin),FC%mat(iin,jin),ET%mat(iin,jin),& ! Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),& ! Ge%mat(iin,jin),Ta_prec%mat(iin,jin),Ts_prec%mat(iin,jin)) ! Ta_prec%mat(iin,jin)=temperatura_oss%mat(iin,jin) ! Ts_prec%mat(iin,jin)=Ts%mat(iin,jin) ! IF(ET%mat(iin,jin)<0)THEN ! ET%mat(iin,jin)=0 ! END IF ! IF (PRESENT(stato_neve)) THEN !IF (stato_neve%mat(iin,jin) /= 0.0) THEN ! ET % mat(iin,jin) = 0.0 ! IF(temperatura_oss%mat(iin,jin)<0.)THEN ! Ts%mat(iin,jin)=temperatura_oss%mat(iin,jin) ! ELSE ! Ts%mat(iin,jin)=0. ! END IF ! Rnetta%mat(iin,jin) = 0.0 ! Xle%mat(iin,jin) = 0.0 ! Ge%mat(iin,jin) = 0.0 ! Hse%mat(iin,jin) = 0.0 ! END IF !END IF !ELSE IF ( snowCovered .OR. iceCovered ) THEN et % mat(i,j) = 0.0 ELSE !compute actual evaporation as a fraction alpha of potential sm = soilMoistureRZ % mat(i,j) wp = wiltingPoint % mat(i,j) fc = fieldCapacity % mat(i,j) alpha = 0.082 * sm + 9.173 * sm**2. - 9.815 * sm**3. IF (alpha > 1.) alpha = 1. actualE = pet % mat(i,j) * alpha !compute actual transpiration as a fraction beta of potential beta = (sm - wp) / (fc - wp) IF (beta > 1.) beta = 1. IF (beta < 0.) beta = 0. actualT = pet % mat(i,j) * beta !compute actual evapotranspiration et % mat(i,j) = vf % mat(i,j) * actualT + (1. - vf % mat(i,j)) * actualE IF ( et % mat(i,j) < 0.) THEN et % mat(i,j) = 0.0 END IF END IF !END IF !compute percolation and capilary rise !------------------------------------------------------------------------ CALL PercolationAndCaprise (id, i, j, rainBalance % mat (i, j), vadose, pet, & soilDepthCell, percolationcellRZ, & percolationcellTZ, runoffcell) !update soil moisture !-------------------- SELECT CASE (infiltrationModel) CASE (SCS_CN, PHILIPEQ, GREEN_AMPT) CALL UpdateSoilMoisture (id, i, j, dtSoilBalance, & soilDepthCell, & percolationcellRZ, & percolationcellTZ, & runoffcell, & QinSoilSub % mat (i,j), & QoutSoilSub % mat (i,j) ) CASE (ROSS_BC) !store values of cumulated variables at previous time step pevapcum = evapcum % mat(i,j) prunoffcum = runoffcum % mat(i,j) pinfilcum = infilcum % mat(i,j) pdrncum = drncum % mat(i,j) !update soil moisture profile and cumulated fluxes in cm CALL SolveRichardsBC (idt = dtSoilBalance, & rain = rainBalance % mat (i, j) , & pet = et % mat (i,j), & p = parameters (cell % mat(i,j) ), & Scell = S (cell % mat(i,j),:), & dxcell = dx (cell % mat(i,j),:), & nsol = 0, & h0 = h0grid % mat(i,j), & nsteps = nsteps % mat(i,j) , & evap = evapcum % mat(i,j), & runoff = runoffcum % mat(i,j), & infil = infilcum % mat(i,j), & drn = drncum % mat(i,j), & watercontent = water ) !compute fluxes of this time step in m/s runoffcell = ( runoffcum % mat(i,j) - prunoffcum ) / 100. / dtSoilBalance IF (runoffcell < 0.) runoffcell = 0. infilt % mat (i,j) = ( infilcum % mat(i,j) - pinfilcum ) / 100. / dtSoilBalance IF (infilt % mat (i,j) < 0.) infilt % mat (i,j) = 0. et % mat (i,j) = ( evapcum % mat(i,j) - pevapcum ) / 100. / dtSoilBalance IF (et % mat (i,j) < 0.) et % mat (i,j) = 0. percolationcellTZ = ( drncum % mat(i,j) - pdrncum ) / 100. / dtSoilBalance IF (percolationcellTZ < 0.) percolationcellTZ = 0. !compute average soil moisture soilMoisture % mat (i,j) = water / soilDepthCell END SELECT END SELECT !update percolation grid percolation % mat (i,j) = percolationcellTZ !update runoff grid runoff % mat (i,j) = runoffcell END DO END DO RETURN END SUBROUTINE SolveSoilBalance !============================================================================== !| Description: ! return .TRUE. if wetland is flooded FUNCTION WetlandIsWet & ! (wetland,row,col,time) & RESULT (wet) IMPLICIT NONE ! function arguments ! Scalar arguments with intent(in): TYPE (grid_integer):: wetland INTEGER, INTENT (IN) :: row, col TYPE (DateTime), INTENT (IN) :: time !Local scalar: LOGICAL :: wet INTEGER :: code INTEGER :: i !------------end of declaration------------------------------------------------ IF (ALLOCATED (wetland % mat) ) THEN code = wetland % mat(row,col) wet = .FALSE. DO i = 1, SIZE(wetCode) IF ( code == wetCode (i) ) THEN IF (time >= wetBegin (i) .AND. time <= wetEnd (i) ) THEN wet = .TRUE. EXIT END IF END IF END DO ELSE wet = .FALSE. END IF END FUNCTION WetlandIsWet !============================================================================== !| Description: ! set initial condition for soilwater balance SUBROUTINE SetInitialCondition & ! (iniDB) IMPLICIT NONE ! Arguments with intent(in): TYPE (IniList), INTENT(IN) :: iniDB !Local declaration: INTEGER (KIND = short) :: i, j REAL (KIND = float) :: scalar !------------end of declaration------------------------------------------------ !mandatory variables ! root-zone soil saturation degree IF (SectionIsPresent('saturation-rz', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'saturation-rz') ) THEN scalar = IniReadReal ('scalar', iniDB, 'saturation-rz') CALL NewGrid (soilSatRZ, mask, scalar) ELSE CALL GridByIni (iniDB, soilSatRZ, section = 'saturation-rz') END IF ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading saturation-rz: ' , & argument = 'section not defined in ini file' ) END IF ! transmission-zone soil saturation degree IF (SectionIsPresent('saturation-tz', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'saturation-tz') ) THEN scalar = IniReadReal ('scalar', iniDB, 'saturation-tz') CALL NewGrid (soilSatTZ, mask, scalar) ELSE CALL GridByIni (iniDB, soilSatTZ, section = 'saturation-tz') END IF ELSE !grid is mandatory: stop the program if not present CALL Catch ('error', 'SoilBalance', & 'error in loading saturation-tz: ' , & argument = 'section not defined in ini file' ) END IF ! allocate mean saturation map CALL NewGrid (soilSat, mask, 0.) !allocate and set soil moisture CALL NewGrid (soilMoisture, mask) CALL NewGrid (soilMoistureRZ, mask) CALL NewGrid (soilMoistureTZ, mask) DO j = 1, mask % jdim DO i = 1, mask % idim SELECT CASE ( balanceId % mat (i,j) ) CASE(LAKE) !lake cells are saturated by definition soilSat % mat(i,j) = 1. soilSatRZ % mat(i,j) = 1. soilSatTZ % mat(i,j) = 1. soilMoisture % mat(i,j) = 1. soilMoistureRZ % mat(i,j) = 1. soilMoistureTZ % mat(i,j) = 1. CASE DEFAULT soilMoistureRZ % mat(i,j) = thetar % mat(i,j) + & soilSatRZ % mat(i,j) * & (thetas % mat(i,j) - & thetar % mat(i,j) ) soilMoistureTZ % mat(i,j) = thetar % mat(i,j) + & soilSatTZ % mat(i,j) * & (thetas % mat(i,j) - & thetar % mat(i,j) ) END SELECT END DO END DO !optional variables: ! interstorm duration IF (SectionIsPresent('interstorm-duration', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'interstorm-duration') ) THEN scalar = IniReadReal ('scalar', iniDB, 'interstorm-duration') CALL NewGrid (interstormDuration, mask, INT(scalar)) ELSE CALL GridByIni (iniDB, interstormDuration, section = 'interstorm-duration') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( interstormDuration, mask, 0 ) END IF ! precipitation status IF (SectionIsPresent('precipitation-status', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'precipitation-status') ) THEN scalar = IniReadReal ('scalar', iniDB, 'precipitation-status') CALL NewGrid (rainFlag, mask, INT(scalar)) ELSE CALL GridByIni (iniDB, rainFlag, section = 'precipitation-status') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( rainFlag, mask, 0 ) END IF ! variables for SCS-CN method IF ( infiltrationModel == SCS_CN ) THEN ! effective soil retention capacity of SCS-CN method [mm] IF (SectionIsPresent('soil-retention', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'soil-retention') ) THEN scalar = IniReadReal ('scalar', iniDB, 'soil-retention') CALL NewGrid (sEff, mask, scalar ) ELSE CALL GridByIni (iniDB, sEff, section = 'soil-retention') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( sEff, mask, 0. ) END IF ! cumulative precipitation IF (SectionIsPresent('cumulative-precipitation', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'cumulative-precipitation') ) THEN scalar = IniReadReal ('scalar', iniDB, 'cumulative-precipitation') CALL NewGrid (cuminf, mask, scalar ) ELSE CALL GridByIni (iniDB, cuminf, section = 'cumulative-precipitation') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( cuminf, mask, 0. ) END IF END IF ! variables for Philip and Green-Ampt methods IF ( infiltrationModel == PHILIPEQ .OR. & infiltrationModel == GREEN_AMPT ) THEN ! cumulative infiltration IF (SectionIsPresent('cumulative-infiltration', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'cumulative-infiltration') ) THEN scalar = IniReadReal ('scalar', iniDB, 'cumulative-infiltration') CALL NewGrid (cuminf, mask, scalar ) ELSE CALL GridByIni (iniDB, cuminf, section = 'cumulative-infiltration') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( cuminf, mask, 0. ) END IF END IF !state variable initialization !IF (SectionIsPresent('initial-saturation', iniDB)) THEN ! ! !cold start ! IF ( KeyIsPresent ('cold', iniDB, section = 'initial-saturation') ) THEN ! !allocate state variables ! CALL NewGrid (sEff, mask, 0.) ! CALL NewGrid (rainFlag, mask, 0) ! CALL NewGrid (interstormDuration, mask, 0) ! ! !initial saturation ! isd = IniReadReal ('cold', iniDB, section = 'initial-saturation') ! ! CALL Catch ('info', 'SoilBalance: ', & ! 'initial degree of saturation: ', & ! argument = ToString(isd)) ! ! !same value for root and transmission zones ! CALL NewGrid (soilSat, mask, isd) ! CALL NewGrid (soilSatRZ, mask, isd) ! CALL NewGrid (soilSatTZ, mask, isd) ! ! ! !allocate and set soil moisture ! CALL NewGrid (soilMoisture, mask) ! CALL NewGrid (soilMoistureRZ, mask) ! CALL NewGrid (soilMoistureTZ, mask) ! ! DO i = 1, mask % idim ! DO j = 1, mask % jdim ! SELECT CASE ( balanceId % mat (i,j) ) ! CASE(LAKE) ! soilSat % mat(i,j) = 1. ! soilSatRZ % mat(i,j) = 1. ! soilSatTZ % mat(i,j) = 1. ! soilMoisture % mat(i,j) = 1. ! CASE DEFAULT ! soilMoisture % mat(i,j) = thetar % mat(i,j) + & ! soilSat % mat(i,j) * & ! (thetas % mat(i,j) - & ! thetar % mat(i,j) ) ! ! !lake cells are saturated by definition ! ! END SELECT ! END DO ! END DO ! ! !same initial soil mositure for root and transmission zones ! soilMoistureRZ = soilMoisture ! soilMoistureTZ = soilMoisture ! ! ELSE !hot start ! !soil moisture ! !TODO HOT START FOR ROOT AND TRANSMISSION ZONES ! CALL GridByIni (iniDB, soilMoisture, section = 'initial-saturation') ! IF ( .NOT. CRSisEqual (mask = mask, grid = soilMoisture, & ! checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance ', & ! 'wrong spatial reference in soil-moisture' ) ! END IF ! ! !compute soil relative saturation ! CALL NewGrid (soilSat, mask) ! ! DO i = 1, mask % idim ! DO j = 1, mask % jdim ! SELECT CASE ( balanceId % mat (i,j) ) ! !lake cells are saturated by definition ! CASE(LAKE) ! soilSat % mat(i,j) = 1. ! soilMoisture % mat(i,j) = thetas % mat(i,j) ! CASE DEFAULT ! soilSat % mat(i,j) = ( soilMoisture % mat(i,j) - & ! thetar % mat(i,j)) / & ! ( thetas % mat(i,j) - & ! thetar % mat(i,j) ) ! END SELECT ! END DO ! END DO ! ! ! ! effective soil retention capacity of SCS-CN method [mm] ! IF (infiltrationModel == SCS_CN ) THEN ! IF (SectionIsPresent('soil-retention', iniDB)) THEN ! CALL GridByIni (iniDB, sEff, section = 'soil-retention') ! IF ( .NOT. CRSisEqual (mask = mask, grid = sEff, & ! checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance', & ! 'wrong spatial reference in soil retention capacity sEff' ) ! END IF ! ELSE ! CALL Catch ('error', 'SoilWaterBalance: ', & ! 'missing soil-retention section in configuration file' ) ! END IF ! END IF ! ! ! ! !cumulative infiltration ! IF (infiltrationModel == PHILIPEQ .OR. & ! infiltrationModel == GREEN_AMPT ) THEN ! IF (SectionIsPresent('cumulative-infiltration', iniDB)) THEN ! CALL GridByIni (iniDB, cuminf, section = 'cumulative-infiltration') ! IF ( .NOT. CRSisEqual (mask = mask, grid = cuminf, & ! checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance', & ! 'wrong spatial reference in cumulative infiltration' ) ! END IF ! ELSE ! CALL Catch ('error', 'SoilWaterBalance: ', & ! 'missing cumulative-infiltration section in configuration file' ) ! END IF ! END IF ! ! ! !precipitation status ! IF (SectionIsPresent('precipitation-status', iniDB)) THEN ! CALL GridByIni (iniDB, rainFlag, section = 'precipitation-status') ! IF ( .NOT. CRSisEqual (mask = mask, grid = rainFlag, & ! checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance', & ! 'wrong spatial reference in precipitation status rainFlag' ) ! END IF ! ELSE ! CALL Catch ('error', 'SoilBalance: ', & ! 'missing precipitation-status section in configuration file' ) ! END IF ! ! ! ! !interstorm duration ! IF (SectionIsPresent('interstorm-duration', iniDB)) THEN ! CALL GridByIni (iniDB, interstormDuration, section = 'interstorm-duration') ! IF ( .NOT. CRSisEqual (mask = mask, grid = interstormDuration, & ! checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance', & ! 'wrong spatial reference in interstorm duration' ) ! END IF ! ELSE ! CALL Catch ('error', 'SoilBalance: ', & ! 'missing interstorm-duration section in configuration file' ) ! END IF ! ! END IF !hot start !ELSE ! CALL Catch ('error', 'SoilBalance: ', & ! 'missing initial-saturation section in configuration file' ) !END IF RETURN END SUBROUTINE SetInitialCondition !============================================================================== !| Description: ! set wetland for soilwater balance SUBROUTINE SetWetland & ! (iniDB, inifile) IMPLICIT NONE ! Arguments with intent(in): TYPE (IniList), INTENT(IN) :: iniDB CHARACTER (LEN = *), INTENT(IN) :: inifile !Local declaration: TYPE (Table) :: wetlandTable INTEGER (KIND = short) :: k !------------end of declaration------------------------------------------------ IF (SectionIsPresent('wetland', iniDB)) THEN CALL GridByIni (iniDB, wetland, section = 'wetland') !read table of flooded period CALL TableNew ( inifile, wetlandTable ) !allocate arrays ALLOCATE ( wetCode ( TableGetNrows(wetlandTable) ) ) ALLOCATE ( wetBegin ( TableGetNrows(wetlandTable) ) ) ALLOCATE ( wetEnd ( TableGetNrows(wetlandTable) ) ) !fill in arrays DO k = 1, TableGetNrows(wetlandTable) wetCode (k) = StringToShort ( wetlandTable % col(1) % row(k) ) wetBegin (k) = wetlandTable % col(2) % row(k) wetEnd (k) = wetlandTable % col(3) % row(k) END DO END IF RETURN END SUBROUTINE SetWetland !============================================================================== !| Description: ! compute percolation and capilalry rise SUBROUTINE PercolationAndCaprise & ! (id, i, j, rain, vadose, pet, soilDepthCell, percolationcellRZ, & percolationcellTZ, runoffcell) IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id INTEGER (KIND = short), INTENT(IN) :: i,j REAL (KIND = float), INTENT(IN) :: rain TYPE (grid_real), INTENT(IN) :: vadose !!vadose zone depth TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration !Arguments with intent(out): REAL (KIND = double), INTENT(OUT) :: soilDepthCell REAL (KIND = float), INTENT(OUT) :: percolationcellRZ REAL (KIND = float), INTENT(OUT) :: percolationcellTZ REAL (KIND = double), INTENT(OUT) :: runoffcell !Local declaration: REAL (KIND = float) :: smAdjustedRZ !! used to prevent soil moisture going to zero REAL (KIND = float) :: smAdjustedTZ !! used to prevent soil moisture going to zero REAL (KIND = float) :: conductivityCellRZ !!partial saturation soil hydraulic !conductivity of current cell [m/s] REAL (KIND = float) :: conductivityCellTZ !!partial saturation soil hydraulic !conductivity of current cell [m/s] !conductivity of current cell [m/s] REAL (KIND = float) :: psiCell !!matric potential of current cell [m] REAL (KIND = float) :: meanHydCond !!mean hydraulic conductivity used to !compute capillary rise [m/s] REAL (KIND = float) :: dsdt !------------end of declaration------------------------------------------------ !calculate unsaturated hydraulic conductivity [m/s] !root zone IF (soilMoistureRZ % mat(i,j) <= thetar % mat(i,j)) THEN smAdjustedRZ = thetar % mat(i,j) + 0.001 ELSE smAdjustedRZ = soilMoistureRZ % mat(i,j) END IF conductivityCellRZ = UnsHydCond (ksat = ksat % mat(i,j), & theta = smAdjustedRZ, & thetas = thetas % mat(i,j), & thetar = thetar % mat(i,j), & psdi = psdi % mat(i,j) ) !transmission zone IF (soilMoistureTZ % mat(i,j) <= thetar % mat(i,j)) THEN smAdjustedTZ = thetar % mat(i,j) + 0.001 ELSE smAdjustedTZ = soilMoistureTZ % mat(i,j) END IF conductivityCellTZ = UnsHydCond (ksat = ksat % mat(i,j), & theta = smAdjustedTZ, & thetas = thetas % mat(i,j), & thetar = thetar % mat(i,j), & psdi = psdi % mat(i,j) ) IF (id == LANDPLAIN ) THEN !set local soil depth !soilDepthCell = MIN (soildepth % mat(i,j), vadose % mat(i,j) ) soilDepthCell = soildepth % mat(i,j) !compute root zone soil suction (m) psiCell = Psi (psic = psic % mat(i,j), theta = smAdjustedRZ, & thetas = thetas % mat(i,j), & thetar = thetar % mat(i,j), & psdi = psdi % mat(i,j) ) !Compute harmonic mean among saturated conductivity !at the interface with groundwater table and unsaturated !conductivity of vadose zone meanHydCond = 2. * conductivityCellRZ * ksat % mat(i,j) / & (conductivityCellRZ + Ksat % mat(i,j)) !compute capillary rise and percolation IF ( vadose % mat (i,j) < soilDepthRZ % mat (i,j) ) THEN ! water table depth lies within the root zone or above ground saturatedByGroundwater = .TRUE. ELSE capRise % mat (i,j) = meanHydCond * psiCell / vadose % mat (i,j) dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * & soilDepthRZ % mat(i,j) / dtSoilBalance percolationcellRZ = MIN ( conductivityCellRZ, dsdt) dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * & soilDepthTZ % mat(i,j) / dtSoilBalance percolationcellTZ = MIN ( conductivityCellTZ, dsdt) * & percolationFactor % mat (i,j) saturatedByGroundwater = .FALSE. END IF ELSE !slope or channel cell soilDepthCell = soilDepth % mat(i,j) capRise % mat(i,j) = 0. dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * & soilDepthRZ % mat(i,j) / dtSoilBalance percolationcellRZ = MIN ( conductivityCellRZ, dsdt) dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * & soilDepthTZ % mat(i,j) / dtSoilBalance percolationcellTZ = MIN ( ksat % mat(i,j), dsdt) * & percolationFactor % mat (i,j) END IF RETURN END SUBROUTINE PercolationAndCaprise !============================================================================== !| Description: ! update soil moisture SUBROUTINE UpdateSoilMoisture & ! (id, i, j, dt, soilDepthCell, percolationcellRZ, percolationcellTZ, & runoffcell, Qin, Qout ) IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id INTEGER (KIND = short), INTENT(IN) :: i,j INTEGER (KIND = short), INTENT(IN) :: dt REAL (KIND = double), INTENT(IN) :: soilDepthCell REAL (KIND = float), INTENT(IN) :: Qin REAL (KIND = float), INTENT(IN) :: Qout !Arguments with intent(inout): REAL (KIND = double), INTENT(INOUT) :: runoffcell REAL (KIND = float), INTENT(INOUT) :: percolationcellRZ REAL (KIND = float), INTENT(INOUT) :: percolationcellTZ !Local declaration: REAL (KIND = double) :: smPrevRZ !! root zone soil moisture of the previous time step REAL (KIND = double) :: smPrevTZ !! transmission zone soil moisture of the previous time step REAL (KIND = float) :: check REAL (KIND = double) :: dsRZ !!root zone soil mositure excess rate (m/s) REAL (KIND = double) :: dsTZ !!transmission zone soil mositure excess rate (m/s) REAL (KIND = float) :: sdRZ !!root zone soil depth (m) REAL (KIND = float) :: sdTZ !!transmission zone soil depth (m) REAL (KIND = double) :: subIn, subOut !!subsurface fluxes REAL (KIND = double) :: area !!cell area (m2) !------------end of declaration------------------------------------------------ IF ( id == LANDPLAIN .AND. saturatedByGroundwater ) THEN ! water table depth lies within the root zone or above ground percolationcellRZ = 0. percolationcellTZ = 0. infilt % mat(i,j) = 0. et % mat(i,j) = pet % mat(i,j) capRise % mat (i,j) = et % mat(i,j) runoffcell = rainBalance % mat (i, j) soilMoistureRZ % mat(i,j) = thetas % mat (i,j) soilMoistureTZ % mat(i,j) = thetas % mat (i,j) soilMoisture % mat(i,j) = thetas % mat (i,j) balanceError % mat (i,j) = 0. RETURN END IF IF ( id == LAKE ) THEN soilMoistureRZ % mat(i,j) = thetas % mat (i,j) soilMoistureTZ % mat(i,j) = thetas % mat (i,j) soilMoisture % mat(i,j) = thetas % mat (i,j) balanceError % mat (i,j) = 0. deltaSoilMoisture % mat (i,j) = 0. RETURN END IF smPrevRZ = soilMoistureRZ % mat(i,j) smPrevTZ = soilMoistureTZ % mat(i,j) sdRZ = soilDepthRZ % mat (i,j) sdTZ = soilDepthTZ % mat (i,j) !IF ( soilDepthCell < soilDepthRZ % mat(i,j) ) THEN ! sdRZ = soilDepthCell ! sdTZ = 0. !ELSE ! sdRZ = soilDepthRZ % mat(i,j) ! sdTZ = MAX (0.1, ( soilDepthCell - sdRZ ) ) !END IF IF ( id == HILLSLOPE .OR. id == CHANNEL ) THEN area = CellArea (soilMoisture, i, j) subIn = Qin / area subOut = Qout / area ELSE IF ( id == LANDPLAIN ) THEN subIn = 0. subOut = 0. END IF IF (sdTZ > 0.) THEN soilMoistureTZ % mat(i,j) = smPrevTZ + & ( percolationcellRZ - & percolationcellTZ + & subIn - subOut ) / & sdTZ * dt IF (soilMoistureTZ % mat(i,j) < thetar % mat(i,j)) THEN soilMoistureTZ % mat(i,j) = thetar % mat(i,j) END IF IF ( soilMoistureTZ % mat(i,j) > thetas % mat(i,j) ) THEN !saturation excess dsTZ = ( soilMoistureTZ % mat(i,j) - thetas % mat(i,j) ) * sdTZ / dt soilMoistureTZ % mat(i,j) = thetas % mat(i,j) ELSE dsTZ = 0. END IF END IF IF (sdRZ > 0.) THEN soilMoistureRZ % mat(i,j) = smPrevRZ + ( infilt % mat(i,j) - percolationcellRZ - & et % mat(i,j) + capRise % mat(i,j) + dsTZ ) / sdRZ * dt END IF !particular cases !case 1: soil is too dry, soil moisture goes below residual value IF (soilMoistureRZ % mat(i,j) < thetar % mat(i,j) ) THEN !compute soil moisture deficit (negative) converted to flux (m/s) dsRZ = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * soilDepthCell / dt !adjust evapotranspiration et % mat(i,j) = et % mat(i,j) + dsRZ !check for negative et IF (et % mat(i,j) < 0. ) et % mat(i,j) = 0. !limit soil moisture soilMoistureRZ % mat(i,j) = thetar % mat(i,j) !case 2: soil is too wet, soil moisture goes above saturated value ELSE IF ( soilMoistureRZ % mat(i,j) > thetas % mat(i,j) ) THEN !compute soil moisture excess converted to flux (m/s) dsRZ = ( soilMoistureRZ % mat(i,j) - thetas % mat(i,j) ) * sdRZ / dt !adjust runoff for saturation excess runoff and infiltration rate runoffcell = runoffcell + dsRZ infilt % mat(i,j) = infilt % mat(i,j) - dsRZ !limit soil moisture to saturation value soilMoistureRZ % mat(i,j) = thetas % mat(i,j) !check if infiltration is negative after correction !adjust capillary rise if present (plain cell) IF (infilt % mat(i,j) < 0.) THEN dsRZ = - infilt % mat(i,j) infilt % mat(i,j) = 0. IF (id == LANDPLAIN ) THEN capRise % mat(i,j) = capRise % mat(i,j) - dsRZ END IF END IF END IF !balance error check = infilt % mat(i,j) + & capRise % mat(i,j) - & percolationcellTZ - & et % mat(i,j) + & subIn - subOut - & (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ / dt - & (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ / dt balanceError % mat (i,j) = check * dt * 1000. !soil mositure variation deltaSoilMoisture % mat (i,j) = (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ + & (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ !update relative saturation and mean soil moisture soilSatRZ % mat(i,j) = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) / & ( thetas % mat(i,j) - thetar % mat(i,j) ) soilSatTZ % mat(i,j) = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) / & ( thetas % mat(i,j) - thetar % mat(i,j) ) soilSat % mat(i,j) = ( soilSatRZ % mat(i,j) * sdRZ + & soilSatTZ % mat(i,j) * sdTZ ) / (sdRZ + sdTZ) IF ( sdRZ + sdTZ > 0. ) THEN soilMoisture % mat(i,j) = ( soilMoistureRZ % mat(i,j) * sdRZ + & soilMoistureTZ % mat(i,j) * sdTZ ) / (sdRZ + sdTZ) ELSE soilMoisture % mat(i,j) = thetas % mat (i,j) END IF RETURN END SUBROUTINE UpdateSoilMoisture END MODULE SoilBalance